Coral Cover trends in the ETP
- Coral cover data
- Generalized Linear Mixed Model design and hypothesis testing
  - Coral cover multiple comparison between years
FIGURE 1. Illustration of the influence of most intense ENSO events on the coral cover in the ETP.
- FIGURE 1a. ENSO and Coral Cover
DHW and coral cover change
Supplementary figures
- FIGURE S2: a) Number of regions surveyed per year, b): Number of countries surveyed per year
- FIGURE S3: Variation in the live coral cover for 1970–2014 per thermal regime

R Code for GCB-19-0810

Title: Coral reef resilience to thermal stress in the Eastern Tropical Pacific
Authors: Mauricio Romero-Torres, Alberto Acosta, Eric A. Treml, Fernando A. Zapata, David A. Paz-García, Ana M. Palacio-Castro & James W. Porter
Version 3.0

rm(list = ls())  # Remove all objects from memory

Coral Cover trends in the ETP

Coral cover data

Select coral cover data and convert variable types

# Select data
  cover <- read.csv("Data/Coral_Cover_ETP.csv", header=TRUE) # Select file: Coral_Cover_ETP.csv
  
# Convert to numeric variable the response variable and convert to nominal variables the factors
  cover$Coral_cover <- as.numeric(cover$Coral_cover) # convert to numeric variable
  cover$Region <- as.factor(cover$Region)               # convert to nominal factors
  cover$Location <- as.factor(cover$Location)
  cover$Site <- as.factor(cover$Site) 
  cover$Country <- as.factor(cover$Country) 
  cover$Lat <- as.numeric(cover$Lat)
  cover$Lon <- as.numeric(cover$Lon)
  cover$Year <- as.numeric((cover$Year)-1970)  # transform years to begin in 0
  summary(cover)

##       ID1                          Reference                    Province  
##  Min.   :  1.0   Fong, 2017             : 33   Cortezian            : 10  
##  1st Qu.:144.2   Glynn, 1983            : 31   Galapagos            : 61  
##  Median :285.5   Calderón-Aguilera, 2006: 26   Tropical_East_Pacific:495  
##  Mean   :285.1   Glynn, 1990            : 24                              
##  3rd Qu.:426.8   Guzman, 2004           : 24                              
##  Max.   :568.0   Cortés, 2010           : 23                              
##                  (Other)                :405                              
##                       Ecoregion                       Region   
##  Nicoya                    :238   Gulf_of_Chiriqui       :133  
##  Panama_Bight              :117   Colombia               : 75  
##  Cortezian                 : 71   Mexican_Central        : 72  
##  Mexican_Tropical_Pacific  : 50   Costa_Rica_Southern    : 56  
##  Eastern_Galapagos_Islands : 31   Nicaragua_Papagayo_Zone: 49  
##  Northern_Galapagos_Islands: 28   Panama_Bight           : 42  
##  (Other)                   : 31   (Other)                :139  
##        Country           Platform             Thermal_regime
##  Panama    :151   Continental:470   Galapagos_Islands: 61   
##  Costa_Rica:143   Oceanic    : 96   Seasonal         : 78   
##  Mexico    :123                     Thermally_Stable :250   
##  Colombia  : 75                     Upwelling        :177   
##  Ecuador   : 61                                             
##  Nicaragua : 11                                             
##  (Other)   :  2                                             
##              Location               Site          Year            Date    
##  Gulf_of_Chiriqui:121   Uva_Island    : 67   Min.   : 0.00   1/1/09 : 55  
##  Gorgona_Island  : 56   La_Azufrada   : 29   1st Qu.:17.00   1/1/06 : 32  
##  Pearl_Islands   : 40   Cano_island   : 27   Median :29.00   1/1/02 : 31  
##  Bahia_Banderas  : 33   Gorgona_Island: 26   Mean   :26.55   1/1/03 : 31  
##  Gulf_of_Papagayo: 30   Saboga_Island : 24   3rd Qu.:36.00   1/1/97 : 26  
##  Los_Cabos       : 28   Cabo_Pulmo    : 17   Max.   :44.00   1/1/98 : 22  
##  (Other)         :258   (Other)       :376                   (Other):369  
##     Period_ENOS   Coral_cover      Algal_cover         Lat        
##  1971-1972: 24   Min.   : 0.000   Min.   : 0.00   Min.   :-1.352  
##  1973-1982: 60   1st Qu.: 8.568   1st Qu.:12.61   1st Qu.: 5.551  
##  1983-1997:169   Median :25.000   Median :35.84   Median : 8.578  
##  1998-2016:313   Mean   :29.185   Mean   :37.57   Mean   : 9.499  
##                  3rd Qu.:45.725   3rd Qu.:54.81   3rd Qu.:10.791  
##                  Max.   :95.720   Max.   :98.00   Max.   :25.796  
##                                   NA's   :390                     
##       Lon         
##  Min.   :-111.23  
##  1st Qu.: -91.82  
##  Median : -83.87  
##  Mean   : -87.92  
##  3rd Qu.: -81.75  
##  Max.   : -77.34  
##

Hierarchically aggregate by Site, Location and Region
Create time-intervals, and
Exclude regions with cero and only one observation of coral cover per time-interval

# Aggregate by site
  aggr.site <- aggregate(Coral_cover ~ Year+Region+Location+Site+Period_ENOS+Thermal_regime,FUN=mean,data=cover)
# Aggregate by location
  aggr.location <- aggregate(Coral_cover~Year+Region+Location+Period_ENOS+Thermal_regime,FUN=mean,data=aggr.site)
# Aggregate by region
  aggr.region <- aggregate(Coral_cover~Year+Region+Period_ENOS+Thermal_regime,FUN=mean,data=aggr.location)

# Create groups for each of three ENSO time-intervals  
  ENSO.73_82 <- subset(aggr.region, Period_ENOS == "1973-1982")
  ENSO.83_97 <- subset(aggr.region, Period_ENOS == "1983-1997")
  ENSO.98_14 <- subset(aggr.region, Period_ENOS == "1998-2016")

# Exclude regions with cero and only one observation of coral cover per time-interval
  ENSO.83_97<- ENSO.83_97[!ENSO.83_97$Region == "Costa_Rica_Central", ]  
  ENSO.83_97<- ENSO.83_97[!ENSO.83_97$Region == "El_Salvador", ]  
  ENSO.83_97<- ENSO.83_97[!ENSO.83_97$Region == "Western_Galapagos_Islands", ]  
  ENSO.83_97<- ENSO.83_97[!ENSO.83_97$Region == "Clipperton", ]  
  ENSO.83_97<- ENSO.83_97[!ENSO.83_97$Region == "Mexican_Southern", ]  
  ENSO.83_97<- ENSO.83_97[!ENSO.83_97$Region == "Northern_Galapagos_Islands", ]  
  ENSO.98_14<- ENSO.98_14[!ENSO.98_14$Region == "Clipperton", ]  
  ENSO.98_14<- ENSO.98_14[!ENSO.98_14$Region == "El_Salvador", ]  
  ENSO.98_14<- ENSO.98_14[!ENSO.98_14$Region == "Western_Galapagos_Islands", ]  
  ENSO.98_14<- ENSO.98_14[!ENSO.98_14$Region == "Revillagigedo_Islands", ]  
  
  ENSO.83_14 <- rbind(ENSO.83_97,ENSO.98_14)   # Second and third time intervals
  ENSO.73_14 <- rbind(ENSO.73_82,ENSO.83_14)   # First, second and third time intervals
  ENSO.73_97 <- rbind(ENSO.73_82,ENSO.83_97)   # First and second and time intervals

Generalized Linear Mixed Model design and hypothesis testing

library(nlme)
# Coral_cover     = response variable
# Region          = repeated measure "subject" or repeated experimental unit (random factor)
# Year            = continuous numeric variable
# Thermal_regime  = fixed factor

Model 1: All years 1970_2014

# All years 1970_2014 
  model1 <- lme(Coral_cover ~ -1 + Thermal_regime * Year, random = ~1|Region, data=aggr.region) 
  summary(model1)

## Linear mixed-effects model fit by REML
##  Data: aggr.region 
##        AIC     BIC    logLik
##   1727.902 1760.58 -853.9508
## 
## Random effects:
##  Formula: ~1 | Region
##         (Intercept) Residual
## StdDev:     12.5915 16.52254
## 
## Fixed effects: Coral_cover ~ -1 + Thermal_regime * Year 
##                                        Value Std.Error  DF   t-value
## Thermal_regimeGalapagos_Islands     21.35928 11.933227  14  1.789900
## Thermal_regimeSeasonal              46.11723 14.871840 181  3.100977
## Thermal_regimeThermally_Stable      17.76232  7.784518 181  2.281749
## Thermal_regimeUpwelling             28.33048  8.794842  14  3.221260
## Year                                -0.11837  0.287082 181 -0.412314
## Thermal_regimeSeasonal:Year         -0.38447  0.509939 181 -0.753959
## Thermal_regimeThermally_Stable:Year  0.34479  0.343646 181  1.003329
## Thermal_regimeUpwelling:Year         0.29697  0.342000 181  0.868347
##                                     p-value
## Thermal_regimeGalapagos_Islands      0.0951
## Thermal_regimeSeasonal               0.0022
## Thermal_regimeThermally_Stable       0.0237
## Thermal_regimeUpwelling              0.0062
## Year                                 0.6806
## Thermal_regimeSeasonal:Year          0.4519
## Thermal_regimeThermally_Stable:Year  0.3170
## Thermal_regimeUpwelling:Year         0.3864
##  Correlation: 
##                                     Th_G_I Thrm_S Th_T_S Thrm_U Year  
## Thermal_regimeSeasonal               0.000                            
## Thermal_regimeThermally_Stable       0.000  0.043                     
## Thermal_regimeUpwelling              0.000  0.000  0.000              
## Year                                -0.681  0.000  0.000  0.000       
## Thermal_regimeSeasonal:Year          0.384 -0.726  0.032  0.000 -0.563
## Thermal_regimeThermally_Stable:Year  0.569  0.009 -0.386  0.000 -0.835
## Thermal_regimeUpwelling:Year         0.572  0.000  0.000 -0.349 -0.839
##                                     Th_S:Y T_T_S:
## Thermal_regimeSeasonal                           
## Thermal_regimeThermally_Stable                   
## Thermal_regimeUpwelling                          
## Year                                             
## Thermal_regimeSeasonal:Year                      
## Thermal_regimeThermally_Stable:Year  0.477       
## Thermal_regimeUpwelling:Year         0.473  0.701
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -2.323212218 -0.685668507 -0.007907488  0.578042210  3.071703318 
## 
## Number of Observations: 202
## Number of Groups: 16

  anova(model1)

##                     numDF denDF   F-value p-value
## Thermal_regime          4   179 14.698982  <.0001
## Year                    1   179  0.717103  0.3982
## Thermal_regime:Year     3   179  1.093597  0.3532

  plot(ranef(model1))    # shows symmetrical scatter effects around zero

  plot(model1)           # plot residuals  vs fitted

  resnorm1 <- resid(model1)
  hist(resnorm1, xlab = "Residuals", main = "") # residuals are normally distributed?

  coef.m1 <- as.data.frame(coef(summary(model1)))    # Coefficients of the model
  
  plot(model1, resid(., type = "p") ~ fitted(.) | Region, abline = 0)

  plot(model1, Coral_cover ~ fitted(.) | Region, abline = c(0,1))

  plot(model1, Coral_cover ~ fitted(.) | Thermal_regime, abline = c(0,1))

Model 2: Time-interval 1 - ENSO.73_82

# Time-interval 1 - ENSO.73_82
    model2 <- lme(Coral_cover ~ -1 + Thermal_regime * Year, random = ~1|Region, data=ENSO.73_82) 
    summary(model2)

## Linear mixed-effects model fit by REML
##  Data: ENSO.73_82 
##        AIC      BIC    logLik
##   164.3022 170.4829 -74.15108
## 
## Random effects:
##  Formula: ~1 | Region
##         (Intercept) Residual
## StdDev:    20.41219 11.23257
## 
## Fixed effects: Coral_cover ~ -1 + Thermal_regime * Year 
##                                         Value Std.Error DF    t-value
## Thermal_regimeGalapagos_Islands      7.102519 22.227571  4  0.3195364
## Thermal_regimeThermally_Stable      23.101084 20.098026  4  1.1494205
## Thermal_regimeUpwelling             -4.938454 19.416560  4 -0.2543424
## Year                                 1.682528  1.727903 13  0.9737397
## Thermal_regimeThermally_Stable:Year -1.381424  2.231710 13 -0.6189979
## Thermal_regimeUpwelling:Year         2.530097  2.488405 13  1.0167545
##                                     p-value
## Thermal_regimeGalapagos_Islands      0.7653
## Thermal_regimeThermally_Stable       0.3144
## Thermal_regimeUpwelling              0.8118
## Year                                 0.3480
## Thermal_regimeThermally_Stable:Year  0.5466
## Thermal_regimeUpwelling:Year         0.3278
##  Correlation: 
##                                     Th_G_I Th_T_S Thrm_U Year   T_T_S:
## Thermal_regimeThermally_Stable       0.000                            
## Thermal_regimeUpwelling              0.000  0.000                     
## Year                                -0.709  0.000  0.000              
## Thermal_regimeThermally_Stable:Year  0.549 -0.419  0.000 -0.774       
## Thermal_regimeUpwelling:Year         0.493  0.000 -0.541 -0.694  0.538
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -1.4650132 -0.4802879 -0.2103293  0.5465022  1.8906401 
## 
## Number of Observations: 22
## Number of Groups: 7

    anova(model2)

##                     numDF denDF  F-value p-value
## Thermal_regime          3     4 3.434033  0.1322
## Year                    1    13 3.583180  0.0808
## Thermal_regime:Year     2    13 1.472433  0.2652

    plot(ranef(model2))

    plot(model2)

    resnorm2 <- resid(model2)
    hist(resnorm2, xlab = "Residuals", main = "")

    coef.m2 <- as.data.frame(coef(summary(model2)))    
    
    plot(model2, resid(., type = "p") ~ fitted(.) | Region, abline = 0)

    plot(model2, Coral_cover ~ fitted(.) | Region, abline = c(0,1))

    plot(model2, Coral_cover ~ fitted(.) | Thermal_regime, abline = c(0,1))

    # Note that this first time interval has low number of regions at each year, even the 
    # GLMM model is adjusted, please warning with interpretation/generalization of results

Model 3: Time-interval 3 - ENSO.83_97

# Time-interval 3 - ENSO.83_97
    model3 <- lme(Coral_cover ~ -1 + Thermal_regime * Year, random = ~1|Region, data=ENSO.83_97) 
    summary(model3)

## Linear mixed-effects model fit by REML
##  Data: ENSO.83_97 
##        AIC      BIC    logLik
##   476.7179 497.1484 -228.3589
## 
## Random effects:
##  Formula: ~1 | Region
##         (Intercept) Residual
## StdDev:    4.741892 9.686501
## 
## Fixed effects: Coral_cover ~ -1 + Thermal_regime * Year 
##                                         Value Std.Error DF   t-value
## Thermal_regimeGalapagos_Islands       1.66897 26.009168  8  0.064169
## Thermal_regimeSeasonal               99.66347 22.231744 50  4.482935
## Thermal_regimeThermally_Stable      -11.52841  9.972345 50 -1.156038
## Thermal_regimeUpwelling             -32.89816 10.298430  8 -3.194483
## Year                                 -0.05590  1.415937 50 -0.039477
## Thermal_regimeSeasonal:Year          -2.98148  1.718591 50 -1.734840
## Thermal_regimeThermally_Stable:Year   1.33736  1.490465 50  0.897275
## Thermal_regimeUpwelling:Year          3.64374  1.483632 50  2.455962
##                                     p-value
## Thermal_regimeGalapagos_Islands      0.9504
## Thermal_regimeSeasonal               0.0000
## Thermal_regimeThermally_Stable       0.2532
## Thermal_regimeUpwelling              0.0127
## Year                                 0.9687
## Thermal_regimeSeasonal:Year          0.0889
## Thermal_regimeThermally_Stable:Year  0.3739
## Thermal_regimeUpwelling:Year         0.0176
##  Correlation: 
##                                     Th_G_I Thrm_S Th_T_S Thrm_U Year  
## Thermal_regimeSeasonal               0.000                            
## Thermal_regimeThermally_Stable       0.000  0.004                     
## Thermal_regimeUpwelling              0.000  0.000  0.000              
## Year                                -0.969  0.000  0.000  0.000       
## Thermal_regimeSeasonal:Year          0.798 -0.557 -0.005  0.000 -0.824
## Thermal_regimeThermally_Stable:Year  0.921 -0.002 -0.294  0.000 -0.950
## Thermal_regimeUpwelling:Year         0.925  0.000  0.000 -0.281 -0.954
##                                     Th_S:Y T_T_S:
## Thermal_regimeSeasonal                           
## Thermal_regimeThermally_Stable                   
## Thermal_regimeUpwelling                          
## Year                                             
## Thermal_regimeSeasonal:Year                      
## Thermal_regimeThermally_Stable:Year  0.785       
## Thermal_regimeUpwelling:Year         0.786  0.907
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.99662705 -0.53542092 -0.06265659  0.58765314  2.99408384 
## 
## Number of Observations: 65
## Number of Groups: 10

    anova(model3)

##                     numDF denDF  F-value p-value
## Thermal_regime          4    48 62.03217  <.0001
## Year                    1    48 39.45391  <.0001
## Thermal_regime:Year     3    48 14.56276  <.0001

    plot(ranef(model3))

    plot(model3)

    resnorm3 <- resid(model3)
    hist(resnorm3, xlab = "Residuals", main = "")

    coef.m3 <- as.data.frame(coef(summary(model3)))    
    
    plot(model3, resid(., type = "p") ~ fitted(.) | Region, abline = 0)

    plot(model3, Coral_cover ~ fitted(.) | Region, abline = c(0,1))

    plot(model3, Coral_cover ~ fitted(.) | Thermal_regime, abline = c(0,1))

Model 4: Time-interval 4 - ENSO.98_14

# Time-interval 4 - ENSO.98_14
    model4 <- lme(Coral_cover ~ -1 + Thermal_regime * Year, random = ~1|Region, data=ENSO.98_14) 
    summary(model4)

## Linear mixed-effects model fit by REML
##  Data: ENSO.98_14 
##        AIC      BIC    logLik
##   866.8707 892.5142 -423.4354
## 
## Random effects:
##  Formula: ~1 | Region
##         (Intercept) Residual
## StdDev:    11.83418 15.29875
## 
## Fixed effects: Coral_cover ~ -1 + Thermal_regime * Year 
##                                        Value Std.Error DF   t-value
## Thermal_regimeGalapagos_Islands     -9.97422  52.62312 10 -0.189541
## Thermal_regimeSeasonal              37.44946  29.75231 87  1.258708
## Thermal_regimeThermally_Stable      -7.42166  25.28540 87 -0.293516
## Thermal_regimeUpwelling             88.18706  20.39094 10  4.324815
## Year                                 0.94979   1.43342 87  0.662603
## Thermal_regimeSeasonal:Year         -1.05637   1.63478 87 -0.646184
## Thermal_regimeThermally_Stable:Year  0.00176   1.59369 87  0.001105
## Thermal_regimeUpwelling:Year        -2.42531   1.52709 87 -1.588187
##                                     p-value
## Thermal_regimeGalapagos_Islands      0.8535
## Thermal_regimeSeasonal               0.2115
## Thermal_regimeThermally_Stable       0.7698
## Thermal_regimeUpwelling              0.0015
## Year                                 0.5093
## Thermal_regimeSeasonal:Year          0.5199
## Thermal_regimeThermally_Stable:Year  0.9991
## Thermal_regimeUpwelling:Year         0.1159
##  Correlation: 
##                                     Th_G_I Thrm_S Th_T_S Thrm_U Year  
## Thermal_regimeSeasonal               0.000                            
## Thermal_regimeThermally_Stable       0.000 -0.014                     
## Thermal_regimeUpwelling              0.000  0.000  0.000              
## Year                                -0.984  0.000  0.000  0.000       
## Thermal_regimeSeasonal:Year          0.863 -0.463 -0.002  0.000 -0.877
## Thermal_regimeThermally_Stable:Year  0.885  0.012 -0.425  0.000 -0.899
## Thermal_regimeUpwelling:Year         0.923  0.000  0.000 -0.327 -0.939
##                                     Th_S:Y T_T_S:
## Thermal_regimeSeasonal                           
## Thermal_regimeThermally_Stable                   
## Thermal_regimeUpwelling                          
## Year                                             
## Thermal_regimeSeasonal:Year                      
## Thermal_regimeThermally_Stable:Year  0.790       
## Thermal_regimeUpwelling:Year         0.823  0.844
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -1.98327607 -0.50110676 -0.02557425  0.30976474  2.71954447 
## 
## Number of Observations: 104
## Number of Groups: 12

    anova(model4)

##                     numDF denDF   F-value p-value
## Thermal_regime          4    85 15.541089  <.0001
## Year                    1    85  1.227329  0.2711
## Thermal_regime:Year     3    85  2.983540  0.0358

    plot(ranef(model4))

    plot(model4)

    resnorm4 <- resid(model4)
    hist(resnorm4, xlab = "Residuals", main = "")

    coef.m4 <- as.data.frame(coef(summary(model4)))    
    
    plot(model4, resid(., type = "p") ~ fitted(.) | Region, abline = 0)

    plot(model4, Coral_cover ~ fitted(.) | Region, abline = c(0,1))

    plot(model4, Coral_cover ~ fitted(.) | Thermal_regime, abline = c(0,1))

# Competing model with autocorrelation structure
    model4b <- lme(Coral_cover ~ -1 + Thermal_regime * Year, random = ~1|Region, 
                   data=ENSO.98_14, correlation=corCompSymm(form=~Year|Region)) 
    
    anova(model4, model4b)  # ~ same results

##         Model df      AIC      BIC    logLik   Test L.Ratio p-value
## model4      1 10 866.8707 892.5142 -423.4354                       
## model4b     2 11 868.8707 897.0786 -423.4354 1 vs 2       0       1

Model 5: Time-interval 5 - ENSO.83_14

# Time-interval 5 - ENSO.83_14  
    model5 <- lme(Coral_cover ~ -1 + Thermal_regime * Year, random = ~1|Region, data=ENSO.83_14) 
    summary(model5)

## Linear mixed-effects model fit by REML
##  Data: ENSO.83_14 
##        AIC      BIC    logLik
##   1423.169 1453.983 -701.5844
## 
## Random effects:
##  Formula: ~1 | Region
##         (Intercept) Residual
## StdDev:    12.53351 15.47206
## 
## Fixed effects: Coral_cover ~ -1 + Thermal_regime * Year 
##                                        Value Std.Error  DF   t-value
## Thermal_regimeGalapagos_Islands     -0.09635 17.626528  11 -0.005466
## Thermal_regimeSeasonal              46.57521 14.105900 151  3.301825
## Thermal_regimeThermally_Stable      -4.10754  9.822261 151 -0.418187
## Thermal_regimeUpwelling             34.46192 10.629603  11  3.242071
## Year                                 0.67715  0.453313 151  1.493774
## Thermal_regimeSeasonal:Year         -1.17692  0.601829 151 -1.955565
## Thermal_regimeThermally_Stable:Year  0.20241  0.521173 151  0.388367
## Thermal_regimeUpwelling:Year        -0.67798  0.519259 151 -1.305672
##                                     p-value
## Thermal_regimeGalapagos_Islands      0.9957
## Thermal_regimeSeasonal               0.0012
## Thermal_regimeThermally_Stable       0.6764
## Thermal_regimeUpwelling              0.0078
## Year                                 0.1373
## Thermal_regimeSeasonal:Year          0.0524
## Thermal_regimeThermally_Stable:Year  0.6983
## Thermal_regimeUpwelling:Year         0.1936
##  Correlation: 
##                                     Th_G_I Thrm_S Th_T_S Thrm_U Year  
## Thermal_regimeSeasonal               0.000                            
## Thermal_regimeThermally_Stable       0.000  0.042                     
## Thermal_regimeUpwelling              0.000  0.000  0.000              
## Year                                -0.834  0.000  0.000  0.000       
## Thermal_regimeSeasonal:Year          0.628 -0.572  0.022  0.000 -0.753
## Thermal_regimeThermally_Stable:Year  0.725  0.010 -0.395  0.000 -0.870
## Thermal_regimeUpwelling:Year         0.728  0.000  0.000 -0.379 -0.873
##                                     Th_S:Y T_T_S:
## Thermal_regimeSeasonal                           
## Thermal_regimeThermally_Stable                   
## Thermal_regimeUpwelling                          
## Year                                             
## Thermal_regimeSeasonal:Year                      
## Thermal_regimeThermally_Stable:Year  0.660       
## Thermal_regimeUpwelling:Year         0.658  0.759
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.52794002 -0.47964749 -0.06607999  0.50350339  3.04770181 
## 
## Number of Observations: 169
## Number of Groups: 13

    anova(model5)

##                     numDF denDF   F-value p-value
## Thermal_regime          4   149 14.198648  <.0001
## Year                    1   149  4.270123  0.0405
## Thermal_regime:Year     3   149  3.799303  0.0116

    plot(ranef(model5))

    plot(model5)

    resnorm5 <- resid(model5)
    hist(resnorm5, xlab = "Residuals", main = "")

    coef.m5 <- as.data.frame(coef(summary(model5))) 
    
    plot(model5, resid(., type = "p") ~ fitted(.) | Region, abline = 0)

    plot(model5, Coral_cover ~ fitted(.) | Region, abline = c(0,1))

    plot(model5, Coral_cover ~ fitted(.) | Thermal_regime, abline = c(0,1))

Coral cover multiple comparison between years

Adjusted p-values for all-pairwise comparisons in a one-way layout
Set up ANOVA model
Set up multiple comparisons object for all-pair comparisons
Multiple comparisons after GLM including interaction terms

library(multcomp)

Year 1982 vs 1983 (codes 12 vs 13)

ENSO.73_97$Year <- as.factor(ENSO.73_97$Year)
comp_years_ENSO.73_97 <- aov(Coral_cover ~ Year + Thermal_regime, data = ENSO.73_97)
m.comp_ENSO.73_97 <- glht(comp_years_ENSO.73_97, linfct = mcp(Year = "Tukey"))
summary(m.comp_ENSO.73_97, test = univariate())

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: aov(formula = Coral_cover ~ Year + Thermal_regime, data = ENSO.73_97)
## 
## Linear Hypotheses:
##               Estimate Std. Error t value Pr(>|t|)    
## 5 - 4 == 0    -3.44519   12.08474  -0.285 0.776577    
## 6 - 4 == 0    11.67379   13.58569   0.859 0.393669    
## 7 - 4 == 0     2.20768   17.29237   0.128 0.898846    
## 8 - 4 == 0   -21.04773   13.64079  -1.543 0.128178    
## 9 - 4 == 0    12.32727   13.64079   0.904 0.369826    
## 10 - 4 == 0    2.61704   12.23990   0.214 0.831431    
## 11 - 4 == 0   10.89268   17.29237   0.630 0.531184    
## 12 - 4 == 0    9.29041   10.82386   0.858 0.394186    
## 13 - 4 == 0  -26.56453   12.25535  -2.168 0.034237 *  
## 14 - 4 == 0  -17.28488   11.32342  -1.526 0.132236    
## 15 - 4 == 0   -6.66044   13.76573  -0.484 0.630289    
## 16 - 4 == 0  -16.73332   11.32342  -1.478 0.144790    
## 17 - 4 == 0   -6.97546   10.29300  -0.678 0.500618    
## 18 - 4 == 0   -5.21764   11.53808  -0.452 0.652776    
## 19 - 4 == 0   -9.53355   10.53984  -0.905 0.369396    
## 20 - 4 == 0  -10.69232   17.29237  -0.618 0.538741    
## 21 - 4 == 0    4.34819   12.74709   0.341 0.734231    
## 22 - 4 == 0    0.52130   10.65659   0.049 0.961150    
## 23 - 4 == 0   -6.10824   11.04278  -0.553 0.582255    
## 24 - 4 == 0    3.00209   11.00761   0.273 0.786014    
## 25 - 4 == 0   25.24915   11.51467   2.193 0.032279 *  
## 26 - 4 == 0    8.32457   11.51467   0.723 0.472565    
## 27 - 4 == 0    0.07463   10.75923   0.007 0.994489    
## 28 - 4 == 0   26.29686   17.32517   1.518 0.134395    
## 6 - 5 == 0    15.11898   13.58569   1.113 0.270282    
## 7 - 5 == 0     5.65287   17.29237   0.327 0.744901    
## 8 - 5 == 0   -17.60254   13.64079  -1.290 0.201934    
## 9 - 5 == 0    15.77246   13.64079   1.156 0.252231    
## 10 - 5 == 0    6.06223   12.23990   0.495 0.622239    
## 11 - 5 == 0   14.33787   17.29237   0.829 0.410365    
## 12 - 5 == 0   12.73559   10.82386   1.177 0.244072    
## 13 - 5 == 0  -23.11934   12.25535  -1.886 0.064156 .  
## 14 - 5 == 0  -13.83970   11.32342  -1.222 0.226485    
## 15 - 5 == 0   -3.21526   13.76573  -0.234 0.816128    
## 16 - 5 == 0  -13.28814   11.32342  -1.174 0.245307    
## 17 - 5 == 0   -3.53027   10.29300  -0.343 0.732834    
## 18 - 5 == 0   -1.77246   11.53808  -0.154 0.878435    
## 19 - 5 == 0   -6.08837   10.53984  -0.578 0.565697    
## 20 - 5 == 0   -7.24713   17.29237  -0.419 0.676669    
## 21 - 5 == 0    7.79338   12.74709   0.611 0.543293    
## 22 - 5 == 0    3.96649   10.65659   0.372 0.711070    
## 23 - 5 == 0   -2.66305   11.04278  -0.241 0.810269    
## 24 - 5 == 0    6.44728   11.00761   0.586 0.560304    
## 25 - 5 == 0   28.69434   11.51467   2.492 0.015533 *  
## 26 - 5 == 0   11.76975   11.51467   1.022 0.310882    
## 27 - 5 == 0    3.51982   10.75923   0.327 0.744717    
## 28 - 5 == 0   29.74205   17.32517   1.717 0.091281 .  
## 7 - 6 == 0    -9.46611   18.35985  -0.516 0.608069    
## 8 - 6 == 0   -32.72152   15.11308  -2.165 0.034435 *  
## 9 - 6 == 0     0.65348   15.11308   0.043 0.965657    
## 10 - 6 == 0   -9.05675   13.81029  -0.656 0.514504    
## 11 - 6 == 0   -0.78111   18.35985  -0.043 0.966209    
## 12 - 6 == 0   -2.38339   12.41250  -0.192 0.848390    
## 13 - 6 == 0  -38.23832   13.92671  -2.746 0.007992 ** 
## 14 - 6 == 0  -28.95867   12.90860  -2.243 0.028641 *  
## 15 - 6 == 0  -18.33423   15.08488  -1.215 0.229053    
## 16 - 6 == 0  -28.40711   12.90860  -2.201 0.031688 *  
## 17 - 6 == 0  -18.64925   12.08019  -1.544 0.127987    
## 18 - 6 == 0  -16.89143   13.27936  -1.272 0.208362    
## 19 - 6 == 0  -21.20734   12.30424  -1.724 0.090021 .  
## 20 - 6 == 0  -22.36611   18.35985  -1.218 0.227994    
## 21 - 6 == 0   -7.32560   14.34376  -0.511 0.611455    
## 22 - 6 == 0  -11.15249   12.45224  -0.896 0.374095    
## 23 - 6 == 0  -17.78203   12.71582  -1.398 0.167222    
## 24 - 6 == 0   -8.67170   12.78086  -0.678 0.500113    
## 25 - 6 == 0   13.57536   13.30638   1.020 0.311792    
## 26 - 6 == 0   -3.34922   13.30638  -0.252 0.802147    
## 27 - 6 == 0  -11.59916   12.56356  -0.923 0.359644    
## 28 - 6 == 0   14.62307   18.62181   0.785 0.435440    
## 8 - 7 == 0   -23.25541   18.23689  -1.275 0.207242    
## 9 - 7 == 0    10.11959   18.23689   0.555 0.581063    
## 10 - 7 == 0    0.40936   17.14223   0.024 0.981029    
## 11 - 7 == 0    8.68500   20.93138   0.415 0.679700    
## 12 - 7 == 0    7.08272   16.43022   0.431 0.667982    
## 13 - 7 == 0  -28.77221   17.29679  -1.663 0.101527    
## 14 - 7 == 0  -19.49257   16.66532  -1.170 0.246846    
## 15 - 7 == 0   -8.86813   18.12710  -0.489 0.626500    
## 16 - 7 == 0  -18.94101   16.66532  -1.137 0.260321    
## 17 - 7 == 0   -9.18314   15.95407  -0.576 0.567075    
## 18 - 7 == 0   -7.42533   16.79651  -0.442 0.660051    
## 19 - 7 == 0  -11.74124   16.16342  -0.726 0.470461    
## 20 - 7 == 0  -12.90000   20.93138  -0.616 0.540067    
## 21 - 7 == 0    2.14051   17.69601   0.121 0.904134    
## 22 - 7 == 0   -1.68638   16.10151  -0.105 0.916942    
## 23 - 7 == 0   -8.31592   16.53313  -0.503 0.616847    
## 24 - 7 == 0    0.79441   16.37629   0.049 0.961474    
## 25 - 7 == 0   23.04147   16.81692   1.370 0.175836    
## 26 - 7 == 0    6.11688   16.81692   0.364 0.717357    
## 27 - 7 == 0   -2.13305   16.22072  -0.132 0.895826    
## 28 - 7 == 0   24.08918   21.30941   1.130 0.262862    
## 9 - 8 == 0    33.37500   14.80072   2.255 0.027861 *  
## 10 - 8 == 0   23.66477   13.52755   1.749 0.085426 .  
## 11 - 8 == 0   31.94041   18.23689   1.751 0.085071 .  
## 12 - 8 == 0   30.33813   12.59792   2.408 0.019176 *  
## 13 - 8 == 0   -5.51680   13.52755  -0.408 0.684882    
## 14 - 8 == 0    3.76284   12.90860   0.291 0.771692    
## 15 - 8 == 0   14.38729   14.93497   0.963 0.339314    
## 16 - 8 == 0    4.31440   12.90860   0.334 0.739393    
## 17 - 8 == 0   14.07227   11.92988   1.180 0.242901    
## 18 - 8 == 0   15.83008   12.90266   1.227 0.224739    
## 19 - 8 == 0   11.51417   12.16281   0.947 0.347668    
## 20 - 8 == 0   10.35541   18.23689   0.568 0.572307    
## 21 - 8 == 0   25.39592   14.02252   1.811 0.075220 .  
## 22 - 8 == 0   21.56903   12.12366   1.779 0.080377 .  
## 23 - 8 == 0   14.93949   12.64762   1.181 0.242259    
## 24 - 8 == 0   24.04982   12.43222   1.934 0.057855 .  
## 25 - 8 == 0   46.29688   12.85667   3.601 0.000651 ***
## 26 - 8 == 0   29.37229   12.85667   2.285 0.025949 *  
## 27 - 8 == 0   21.12236   12.22384   1.728 0.089226 .  
## 28 - 8 == 0   47.34459   18.23689   2.596 0.011882 *  
## 10 - 9 == 0   -9.71023   13.52755  -0.718 0.475706    
## 11 - 9 == 0   -1.43459   18.23689  -0.079 0.937566    
## 12 - 9 == 0   -3.03687   12.59792  -0.241 0.810344    
## 13 - 9 == 0  -38.89180   13.52755  -2.875 0.005611 ** 
## 14 - 9 == 0  -29.61216   12.90860  -2.294 0.025369 *  
## 15 - 9 == 0  -18.98771   14.93497  -1.271 0.208590    
## 16 - 9 == 0  -29.06060   12.90860  -2.251 0.028108 *  
## 17 - 9 == 0  -19.30273   11.92988  -1.618 0.110993    
## 18 - 9 == 0  -17.54492   12.90266  -1.360 0.179070    
## 19 - 9 == 0  -21.86083   12.16281  -1.797 0.077398 .  
## 20 - 9 == 0  -23.01959   18.23689  -1.262 0.211824    
## 21 - 9 == 0   -7.97908   14.02252  -0.569 0.571503    
## 22 - 9 == 0  -11.80597   12.12366  -0.974 0.334134    
## 23 - 9 == 0  -18.43551   12.64762  -1.458 0.150246    
## 24 - 9 == 0   -9.32518   12.43222  -0.750 0.456186    
## 25 - 9 == 0   12.92188   12.85667   1.005 0.318967    
## 26 - 9 == 0   -4.00271   12.85667  -0.311 0.756645    
## 27 - 9 == 0  -12.25264   12.22384  -1.002 0.320265    
## 28 - 9 == 0   13.96959   18.23689   0.766 0.446726    
## 11 - 10 == 0   8.27564   17.14223   0.483 0.631050    
## 12 - 10 == 0   6.67336   11.04000   0.604 0.547848    
## 13 - 10 == 0 -29.18157   12.15792  -2.400 0.019560 *  
## 14 - 10 == 0 -19.90193   11.39114  -1.747 0.085816 .  
## 15 - 10 == 0  -9.27748   13.57665  -0.683 0.497066    
## 16 - 10 == 0 -19.35036   11.39114  -1.699 0.094640 .  
## 17 - 10 == 0  -9.59250   10.28678  -0.933 0.354875    
## 18 - 10 == 0  -7.83468   11.45073  -0.684 0.496522    
## 19 - 10 == 0 -12.15059   10.57357  -1.149 0.255132    
## 20 - 10 == 0 -13.30936   17.14223  -0.776 0.440610    
## 21 - 10 == 0   1.73115   12.71073   0.136 0.892130    
## 22 - 10 == 0  -2.09574   10.51193  -0.199 0.842661    
## 23 - 10 == 0  -8.72528   11.12874  -0.784 0.436158    
## 24 - 10 == 0   0.38505   10.88700   0.035 0.971906    
## 25 - 10 == 0  22.63211   11.42620   1.981 0.052291 .  
## 26 - 10 == 0   5.70753   11.42620   0.500 0.619276    
## 27 - 10 == 0  -2.54241   10.64946  -0.239 0.812138    
## 28 - 10 == 0  23.67982   17.29679   1.369 0.176179    
## 12 - 11 == 0  -1.60228   16.43022  -0.098 0.922644    
## 13 - 11 == 0 -37.45721   17.29679  -2.166 0.034399 *  
## 14 - 11 == 0 -28.17757   16.66532  -1.691 0.096154 .  
## 15 - 11 == 0 -17.55312   18.12710  -0.968 0.336830    
## 16 - 11 == 0 -27.62601   16.66532  -1.658 0.102687    
## 17 - 11 == 0 -17.86814   15.95407  -1.120 0.267263    
## 18 - 11 == 0 -16.11033   16.79651  -0.959 0.341399    
## 19 - 11 == 0 -20.42624   16.16342  -1.264 0.211297    
## 20 - 11 == 0 -21.58500   20.93138  -1.031 0.306643    
## 21 - 11 == 0  -6.54449   17.69601  -0.370 0.712834    
## 22 - 11 == 0 -10.37138   16.10151  -0.644 0.521991    
## 23 - 11 == 0 -17.00092   16.53313  -1.028 0.308009    
## 24 - 11 == 0  -7.89059   16.37629  -0.482 0.631709    
## 25 - 11 == 0  14.35647   16.81692   0.854 0.396729    
## 26 - 11 == 0  -2.56812   16.81692  -0.153 0.879148    
## 27 - 11 == 0 -10.81805   16.22072  -0.667 0.507418    
## 28 - 11 == 0  15.40418   21.30941   0.723 0.472609    
## 13 - 12 == 0 -35.85493   11.10858  -3.228 0.002039 ** 
## 14 - 12 == 0 -26.57529    9.96432  -2.667 0.009861 ** 
## 15 - 12 == 0 -15.95085   12.66577  -1.259 0.212857    
## 16 - 12 == 0 -26.02373    9.96432  -2.612 0.011408 *  
## 17 - 12 == 0 -16.26586    8.81152  -1.846 0.069914 .  
## 18 - 12 == 0 -14.50805   10.30126  -1.408 0.164268    
## 19 - 12 == 0 -18.82396    9.10569  -2.067 0.043104 *  
## 20 - 12 == 0 -19.98272   16.43022  -1.216 0.228746    
## 21 - 12 == 0  -4.94221   11.64016  -0.425 0.672685    
## 22 - 12 == 0  -8.76910    9.26628  -0.946 0.347832    
## 23 - 12 == 0 -15.39864    9.67190  -1.592 0.116706    
## 24 - 12 == 0  -6.28831    9.68212  -0.649 0.518551    
## 25 - 12 == 0  15.95875   10.29949   1.549 0.126617    
## 26 - 12 == 0  -0.96584   10.29949  -0.094 0.925606    
## 27 - 12 == 0  -9.21577    9.39664  -0.981 0.330720    
## 28 - 12 == 0  17.00646   16.56830   1.026 0.308872    
## 14 - 13 == 0   9.27964   11.46179   0.810 0.421414    
## 15 - 13 == 0  19.90409   13.77128   1.445 0.153655    
## 16 - 13 == 0   9.83121   11.46179   0.858 0.394508    
## 17 - 13 == 0  19.58907   10.32922   1.896 0.062797 .  
## 18 - 13 == 0  21.34689   11.38881   1.874 0.065832 .  
## 19 - 13 == 0  17.03098   10.57989   1.610 0.112791    
## 20 - 13 == 0  15.87221   17.29679   0.918 0.362544    
## 21 - 13 == 0  30.91272   12.63195   2.447 0.017394 *  
## 22 - 13 == 0  27.08583   10.55141   2.567 0.012813 *  
## 23 - 13 == 0  20.45629   11.13294   1.837 0.071180 .  
## 24 - 13 == 0  29.56662   10.88394   2.717 0.008643 ** 
## 25 - 13 == 0  51.81368   11.30914   4.582 2.44e-05 ***
## 26 - 13 == 0  34.88910   11.30914   3.085 0.003096 ** 
## 27 - 13 == 0  26.63916   10.64424   2.503 0.015116 *  
## 28 - 13 == 0  52.86139   17.14223   3.084 0.003108 ** 
## 15 - 14 == 0  10.62444   12.96928   0.819 0.415969    
## 16 - 14 == 0   0.55156   10.46569   0.053 0.958147    
## 17 - 14 == 0  10.30942    9.32873   1.105 0.273591    
## 18 - 14 == 0  12.06724   10.69153   1.129 0.263604    
## 19 - 14 == 0   7.75133    9.62096   0.806 0.423668    
## 20 - 14 == 0   6.59257   16.66532   0.396 0.693837    
## 21 - 14 == 0  21.63308   12.00358   1.802 0.076620 .  
## 22 - 14 == 0  17.80619    9.69120   1.837 0.071195 .  
## 23 - 14 == 0  11.17665   10.18459   1.097 0.276922    
## 24 - 14 == 0  20.28698   10.09227   2.010 0.048994 *  
## 25 - 14 == 0  42.53404   10.68062   3.982 0.000190 ***
## 26 - 14 == 0  25.60945   10.68062   2.398 0.019679 *  
## 27 - 14 == 0  17.35952    9.82504   1.767 0.082423 .  
## 28 - 14 == 0  43.58175   16.81001   2.593 0.011990 *  
## 16 - 15 == 0 -10.07288   12.96928  -0.777 0.440455    
## 17 - 15 == 0  -0.31502   12.04166  -0.026 0.979217    
## 18 - 15 == 0   1.44280   13.13743   0.110 0.912922    
## 19 - 15 == 0  -2.87311   12.31769  -0.233 0.816374    
## 20 - 15 == 0  -4.03188   18.12710  -0.222 0.824753    
## 21 - 15 == 0  11.00863   14.26948   0.771 0.443501    
## 22 - 15 == 0   7.18174   12.23634   0.587 0.559497    
## 23 - 15 == 0   0.55221   12.79898   0.043 0.965732    
## 24 - 15 == 0   9.66253   12.59573   0.767 0.446065    
## 25 - 15 == 0  31.90959   13.16352   2.424 0.018430 *  
## 26 - 15 == 0  14.98501   13.16352   1.138 0.259566    
## 27 - 15 == 0   6.73508   12.39278   0.543 0.588856    
## 28 - 15 == 0  32.95730   18.56234   1.775 0.080974 .  
## 17 - 16 == 0   9.75786    9.32873   1.046 0.299827    
## 18 - 16 == 0  11.51568   10.69153   1.077 0.285827    
## 19 - 16 == 0   7.19977    9.62096   0.748 0.457227    
## 20 - 16 == 0   6.04101   16.66532   0.362 0.718281    
## 21 - 16 == 0  21.08151   12.00358   1.756 0.084232 .  
## 22 - 16 == 0  17.25462    9.69120   1.780 0.080153 .  
## 23 - 16 == 0  10.62509   10.18459   1.043 0.301088    
## 24 - 16 == 0  19.73541   10.09227   1.955 0.055266 .  
## 25 - 16 == 0  41.98247   10.68062   3.931 0.000225 ***
## 26 - 16 == 0  25.05789   10.68062   2.346 0.022350 *  
## 27 - 16 == 0  16.80796    9.82504   1.711 0.092385 .  
## 28 - 16 == 0  43.03019   16.81001   2.560 0.013055 *  
## 18 - 17 == 0   1.75782    9.36462   0.188 0.851749    
## 19 - 17 == 0  -2.55809    8.24140  -0.310 0.757354    
## 20 - 17 == 0  -3.71686   15.95407  -0.233 0.816589    
## 21 - 17 == 0  11.32365   10.67732   1.061 0.293225    
## 22 - 17 == 0   7.49676    8.27310   0.906 0.368536    
## 23 - 17 == 0   0.86722    8.79756   0.099 0.921809    
## 24 - 17 == 0   9.97755    8.71342   1.145 0.256800    
## 25 - 17 == 0  32.22461    9.44515   3.412 0.001171 ** 
## 26 - 17 == 0  15.30003    9.44515   1.620 0.110591    
## 27 - 17 == 0   7.05009    8.33792   0.846 0.401222    
## 28 - 17 == 0  33.27232   16.03612   2.075 0.042372 *  
## 19 - 18 == 0  -4.31591    9.62361  -0.448 0.655456    
## 20 - 18 == 0  -5.47467   16.79651  -0.326 0.745622    
## 21 - 18 == 0   9.56583   11.53751   0.829 0.410386    
## 22 - 18 == 0   5.73894    9.60194   0.598 0.552337    
## 23 - 18 == 0  -0.89059   10.05120  -0.089 0.929695    
## 24 - 18 == 0   8.21973    9.94532   0.826 0.411855    
## 25 - 18 == 0  30.46679   10.56650   2.883 0.005483 ** 
## 26 - 18 == 0  13.54221   10.56650   1.282 0.204991    
## 27 - 18 == 0   5.29228    9.57737   0.553 0.582637    
## 28 - 18 == 0  31.51450   16.66974   1.891 0.063602 .  
## 20 - 19 == 0  -1.15876   16.16342  -0.072 0.943091    
## 21 - 19 == 0  13.88175   10.86908   1.277 0.206543    
## 22 - 19 == 0  10.05485    8.60023   1.169 0.247050    
## 23 - 19 == 0   3.42532    9.05878   0.378 0.706697    
## 24 - 19 == 0  12.53564    9.01126   1.391 0.169416    
## 25 - 19 == 0  34.78271    9.70911   3.582 0.000690 ***
## 26 - 19 == 0  17.85812    9.70911   1.839 0.070902 .  
## 27 - 19 == 0   9.60819    8.63586   1.113 0.270396    
## 28 - 19 == 0  35.83042   16.17582   2.215 0.030629 *  
## 21 - 20 == 0  15.04051   17.69601   0.850 0.398796    
## 22 - 20 == 0  11.21362   16.10151   0.696 0.488893    
## 23 - 20 == 0   4.58408   16.53313   0.277 0.782545    
## 24 - 20 == 0  13.69441   16.37629   0.836 0.406398    
## 25 - 20 == 0  35.94147   16.81692   2.137 0.036736 *  
## 26 - 20 == 0  19.01688   16.81692   1.131 0.262707    
## 27 - 20 == 0  10.76695   16.22072   0.664 0.509418    
## 28 - 20 == 0  36.98918   21.30941   1.736 0.087817 .  
## 22 - 21 == 0  -3.82689   10.87482  -0.352 0.726164    
## 23 - 21 == 0 -10.45643   10.97923  -0.952 0.344788    
## 24 - 21 == 0  -1.34610   11.13470  -0.121 0.904187    
## 25 - 21 == 0  20.90096   11.88976   1.758 0.083952 .  
## 26 - 21 == 0   3.97638   11.88976   0.334 0.739236    
## 27 - 21 == 0  -4.27356   10.64946  -0.401 0.689654    
## 28 - 21 == 0  21.94867   17.52595   1.252 0.215383    
## 23 - 22 == 0  -6.62954    9.14341  -0.725 0.471280    
## 24 - 22 == 0   2.48079    8.97050   0.277 0.783092    
## 25 - 22 == 0  24.72785    9.68708   2.553 0.013298 *  
## 26 - 22 == 0   7.80327    9.68708   0.806 0.423747    
## 27 - 22 == 0  -0.44667    8.60178  -0.052 0.958762    
## 28 - 22 == 0  25.77556   16.17879   1.593 0.116466    
## 24 - 23 == 0   9.11033    9.50529   0.958 0.341748    
## 25 - 23 == 0  31.35739   10.30850   3.042 0.003505 ** 
## 26 - 23 == 0  14.43280   10.30850   1.400 0.166724    
## 27 - 23 == 0   6.18287    9.04464   0.684 0.496907    
## 28 - 23 == 0  32.40510   16.54163   1.959 0.054844 .  
## 25 - 24 == 0  22.24706   10.03707   2.216 0.030526 *  
## 26 - 24 == 0   5.32248   10.03707   0.530 0.597907    
## 27 - 24 == 0  -2.92746    8.99241  -0.326 0.745919    
## 28 - 24 == 0  23.29477   16.37018   1.423 0.160002    
## 26 - 25 == 0 -16.92458   10.46569  -1.617 0.111181    
## 27 - 25 == 0 -25.17452    9.77607  -2.575 0.012548 *  
## 28 - 25 == 0   1.04771   16.57783   0.063 0.949821    
## 27 - 26 == 0  -8.24993    9.77607  -0.844 0.402139    
## 28 - 26 == 0  17.97229   16.57783   1.084 0.282724    
## 28 - 27 == 0  26.22223   16.21043   1.618 0.111080    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Univariate p values reported)

Year 1983 vs 1997 (codes 13 vs 27), and 1997 vs 1998 (codes 27 vs 28)

ENSO.83_97$Year <- as.factor(ENSO.83_97$Year)
comp_years_ENSO.83_97 <- aov(Coral_cover ~ Year + Thermal_regime, data = ENSO.83_97)
m.comp_ENSO.83_97 <- glht(comp_years_ENSO.83_97, linfct = mcp(Year = "Tukey"))
summary(m.comp_ENSO.83_97, test = univariate())

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: aov(formula = Coral_cover ~ Year + Thermal_regime, data = ENSO.83_97)
## 
## Linear Hypotheses:
##                Estimate Std. Error t value Pr(>|t|)    
## 14 - 13 == 0  12.351578  10.375355   1.190 0.239968    
## 15 - 13 == 0  21.853429  12.451269   1.755 0.085898 .  
## 16 - 13 == 0  12.903140  10.375355   1.244 0.219935    
## 17 - 13 == 0  21.659577   9.319358   2.324 0.024588 *  
## 18 - 13 == 0  21.654675  10.234189   2.116 0.039796 *  
## 19 - 13 == 0  19.121678   9.546497   2.003 0.051090 .  
## 20 - 13 == 0  17.821554  15.603472   1.142 0.259301    
## 21 - 13 == 0  31.408598  11.354258   2.766 0.008138 ** 
## 22 - 13 == 0  27.940801   9.494612   2.943 0.005082 ** 
## 23 - 13 == 0  23.016432  10.063858   2.287 0.026843 *  
## 24 - 13 == 0  30.202718   9.785636   3.086 0.003426 ** 
## 25 - 13 == 0  51.570014  10.163303   5.074 6.85e-06 ***
## 26 - 13 == 0  34.645431  10.163303   3.409 0.001366 ** 
## 27 - 13 == 0  27.374437   9.570674   2.860 0.006346 ** 
## 28 - 13 == 0  51.886723  15.418676   3.365 0.001551 ** 
## 15 - 14 == 0   9.501851  11.732643   0.810 0.422188    
## 16 - 14 == 0   0.551562   9.403948   0.059 0.953483    
## 17 - 14 == 0   9.307999   8.393339   1.109 0.273206    
## 18 - 14 == 0   9.303098   9.673243   0.962 0.341211    
## 19 - 14 == 0   6.770100   8.654255   0.782 0.438055    
## 20 - 14 == 0   5.469976  15.036262   0.364 0.717685    
## 21 - 14 == 0  19.057020  10.836676   1.759 0.085303 .  
## 22 - 14 == 0  15.589224   8.769831   1.778 0.082081 .  
## 23 - 14 == 0  10.664854   9.156204   1.165 0.250118    
## 24 - 14 == 0  17.851140   9.129581   1.955 0.056638 .  
## 25 - 14 == 0  39.218436   9.691117   4.047 0.000197 ***
## 26 - 14 == 0  22.293853   9.691117   2.300 0.026008 *  
## 27 - 14 == 0  15.022859   8.885050   1.691 0.097641 .  
## 28 - 14 == 0  39.535145  15.199628   2.601 0.012458 *  
## 16 - 15 == 0  -8.950288  11.732643  -0.763 0.449445    
## 17 - 15 == 0  -0.193851  10.880330  -0.018 0.985862    
## 18 - 15 == 0  -0.198753  11.878661  -0.017 0.986723    
## 19 - 15 == 0  -2.731750  11.148223  -0.245 0.807516    
## 20 - 15 == 0  -4.031875  16.288116  -0.248 0.805596    
## 21 - 15 == 0   9.555169  12.906122   0.740 0.462844    
## 22 - 15 == 0   6.087373  11.030341   0.552 0.583705    
## 23 - 15 == 0   1.163004  11.599831   0.100 0.920573    
## 24 - 15 == 0   8.349289  11.367353   0.734 0.466374    
## 25 - 15 == 0  29.716586  11.930022   2.491 0.016413 *  
## 26 - 15 == 0  12.792002  11.930022   1.072 0.289200    
## 27 - 15 == 0   5.521009  11.187549   0.493 0.624010    
## 28 - 15 == 0  30.033295  16.807767   1.787 0.080549 .  
## 17 - 16 == 0   8.756437   8.393339   1.043 0.302280    
## 18 - 16 == 0   8.751535   9.673243   0.905 0.370330    
## 19 - 16 == 0   6.218538   8.654255   0.719 0.476051    
## 20 - 16 == 0   4.918413  15.036262   0.327 0.745074    
## 21 - 16 == 0  18.505457  10.836676   1.708 0.094439 .  
## 22 - 16 == 0  15.037661   8.769831   1.715 0.093130 .  
## 23 - 16 == 0  10.113292   9.156204   1.105 0.275108    
## 24 - 16 == 0  17.299577   9.129581   1.895 0.064402 .  
## 25 - 16 == 0  38.666874   9.691117   3.990 0.000235 ***
## 26 - 16 == 0  21.742290   9.691117   2.244 0.029720 *  
## 27 - 16 == 0  14.471297   8.885050   1.629 0.110202    
## 28 - 16 == 0  38.983583  15.199628   2.565 0.013653 *  
## 18 - 17 == 0  -0.004902   8.444989  -0.001 0.999539    
## 19 - 17 == 0  -2.537899   7.407768  -0.343 0.733460    
## 20 - 17 == 0  -3.838024  14.381093  -0.267 0.790755    
## 21 - 17 == 0   9.749021   9.617622   1.014 0.316046    
## 22 - 17 == 0   6.281224   7.459867   0.842 0.404143    
## 23 - 17 == 0   1.356855   7.915355   0.171 0.864645    
## 24 - 17 == 0   8.543141   7.854396   1.088 0.282399    
## 25 - 17 == 0  29.910437   8.539886   3.502 0.001038 ** 
## 26 - 17 == 0  12.985854   8.539886   1.521 0.135202    
## 27 - 17 == 0   5.714860   7.513916   0.761 0.450795    
## 28 - 17 == 0  30.227146  14.473044   2.089 0.042314 *  
## 19 - 18 == 0  -2.532997   8.679641  -0.292 0.771727    
## 20 - 18 == 0  -3.833122  15.150473  -0.253 0.801393    
## 21 - 18 == 0   9.753922  10.369343   0.941 0.351798    
## 22 - 18 == 0   6.286126   8.639110   0.728 0.470524    
## 23 - 18 == 0   1.361757   9.084485   0.150 0.881499    
## 24 - 18 == 0   8.548042   8.940295   0.956 0.344007    
## 25 - 18 == 0  29.915339   9.498092   3.150 0.002871 ** 
## 26 - 18 == 0  12.990755   9.498092   1.368 0.178044    
## 27 - 18 == 0   5.719762   8.609504   0.664 0.509780    
## 28 - 18 == 0  30.232048  14.998026   2.016 0.049691 *  
## 20 - 19 == 0  -1.300125  14.584825  -0.089 0.929356    
## 21 - 19 == 0  12.286919   9.787894   1.255 0.215703    
## 22 - 19 == 0   8.819123   7.765025   1.136 0.261944    
## 23 - 19 == 0   3.894754   8.143592   0.478 0.634730    
## 24 - 19 == 0  11.081039   8.129358   1.363 0.179490    
## 25 - 19 == 0  32.448336   8.775695   3.698 0.000579 ***
## 26 - 19 == 0  15.523753   8.775695   1.769 0.083532 .  
## 27 - 19 == 0   8.252759   7.788304   1.060 0.294843    
## 28 - 19 == 0  32.765045  14.592673   2.245 0.029597 *  
## 21 - 20 == 0  13.587044  15.968790   0.851 0.399260    
## 22 - 20 == 0  10.119248  14.494919   0.698 0.488614    
## 23 - 20 == 0   5.194879  14.932861   0.348 0.729515    
## 24 - 20 == 0  12.381164  14.752999   0.839 0.405680    
## 25 - 20 == 0  33.748461  15.190776   2.222 0.031268 *  
## 26 - 20 == 0  16.823877  15.190776   1.108 0.273833    
## 27 - 20 == 0   9.552884  14.614907   0.654 0.516598    
## 28 - 20 == 0  34.065170  19.259680   1.769 0.083569 .  
## 22 - 21 == 0  -3.467796   9.790236  -0.354 0.724800    
## 23 - 21 == 0  -8.392165   9.902026  -0.848 0.401094    
## 24 - 21 == 0  -1.205880  10.014648  -0.120 0.904681    
## 25 - 21 == 0  20.161416  10.688016   1.886 0.065569 .  
## 26 - 21 == 0   3.236833  10.688016   0.303 0.763372    
## 27 - 21 == 0  -4.034160   9.577556  -0.421 0.675564    
## 28 - 21 == 0  20.478126  15.763214   1.299 0.200381    
## 23 - 22 == 0  -4.924369   8.277348  -0.595 0.554812    
## 24 - 22 == 0   2.261916   8.062383   0.281 0.780312    
## 25 - 22 == 0  23.629213   8.730875   2.706 0.009511 ** 
## 26 - 22 == 0   6.704629   8.730875   0.768 0.446458    
## 27 - 22 == 0  -0.566364   7.731752  -0.073 0.941923    
## 28 - 22 == 0  23.945922  14.586686   1.642 0.107487    
## 24 - 23 == 0   7.186285   8.596189   0.836 0.407485    
## 25 - 23 == 0  28.553582   9.335014   3.059 0.003699 ** 
## 26 - 23 == 0  11.628999   9.335014   1.246 0.219168    
## 27 - 23 == 0   4.358005   8.178387   0.533 0.596690    
## 28 - 23 == 0  28.870291  14.934120   1.933 0.059383 .  
## 25 - 24 == 0  21.367296   9.033108   2.365 0.022280 *  
## 26 - 24 == 0   4.442713   9.033108   0.492 0.625182    
## 27 - 24 == 0  -2.828280   8.080350  -0.350 0.727922    
## 28 - 24 == 0  21.684006  14.744958   1.471 0.148206    
## 26 - 25 == 0 -16.924583   9.403948  -1.800 0.078464 .  
## 27 - 25 == 0 -24.195577   8.799034  -2.750 0.008496 ** 
## 28 - 25 == 0   0.316709  14.905050   0.021 0.983139    
## 27 - 26 == 0  -7.270994   8.799034  -0.826 0.412877    
## 28 - 26 == 0  17.241292  14.905050   1.157 0.253350    
## 28 - 27 == 0  24.512286  14.601377   1.679 0.099979 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Univariate p values reported)

FIGURE 1. Illustration of the influence of most intense ENSO events on the coral cover in the ETP.

# Load packages 
library(dplyr)
library(RColorBrewer)
library(plotrix)
library(readr)
library(reshape)
library(data.table)

FIGURE 1a. ENSO and Coral Cover

# Plot ENOS 1971 - 2014
    ENOS<-read.table(file="http://www.cpc.ncep.noaa.gov/data/indices/ersst5.nino.mth.81-10.ascii",header=TRUE)
    ENOS1 <-ENOS[241:780, ]   # subset from 1982-2014
    Year <- ENOS1$YR
    Anom3 <- ENOS1$ANOM.3
    YR_ANOM3 <- data.frame(cbind(Year,Anom3))  # Create a data frame with two columns
    YR_ANOM3$Year <- as.numeric(YR_ANOM3$Year)
    x <- YR_ANOM3                              # Funcion My Plot to graph ENSO
    Year <- seq(
       as.Date(sprintf("%s/1/1",min(x$Year))), 
       as.Date(sprintf("%s/12/1",max(x$Year))), 
       "months")
   
     myPlot = function(y, x, 
                      colors=c("red","blue"), 
                      xlab=deparse(substitute(x)),
                      ylab=deparse(substitute(y)), ...)
    {
       n = length(x)
       y.pos = y.neg = y
       y.pos[y<0] = 0
       y.neg[y>0] = 0
       plot(y~x, type="n", xlab=xlab, ylab=ylab, ...)   # xaxt="n", Apague eje x
       #grid()
       polygon(c(x[1],x,x[n]),c(0,y.pos,0),col=colors[1],border=NA)
       polygon(c(x[1],x,x[n]),c(0,y.neg,0),col=colors[2],border=NA)
    } 

# Boxplot coral cover all years 
    cover <- read.csv("Data/Coral_Cover_ETP.csv", header=TRUE)      # Select file: Coral_Cover_ETP.csv
    cover$Coral_cover <- as.numeric(cover$Coral_cover)  # convert to numeric variable
    cover$Region <- as.factor(cover$Region)                # convert to nominal factors
    cover$Location <- as.factor(cover$Location)
    cover$Site <- as.factor(cover$Site) 
    cover$Year <- as.numeric(cover$Year)
    cover$Country <- factor(cover$Country) 

# Hierarchically aggregate by Site, Location and Region.
  # Agregate by site
    aggr.site <- aggregate(Coral_cover ~ Year+Region+Location+Site+Period_ENOS+Thermal_regime,FUN=mean,data=cover)
  # Agregate by location
    aggr.location <- aggregate(Coral_cover~Year+Region+Location+Period_ENOS+Thermal_regime,FUN=mean,data=aggr.site)
  # Aggregate by region
    aggr.region <- aggregate(Coral_cover~Year+Region+Period_ENOS+Thermal_regime,FUN=mean,data=aggr.location)

# Plot 
  # dev.off() # reset parameters for plot
  op <- par(mfcol=c(2,1), cex.main = 1.5, mar = c(0,5,0,2), oma=c(6,4,4,2), cex.lab = 1.5, cex.axis = 1,  las = 1)  # bty = "n",
  myPlot(x$A,Year,ylab="",xlab="",axes = FALSE)  
  axis(side=2, at=c(-2,-1,0,1,2), lwd = 0.5)
  box(lwd = 0.5)
  abline(h=0, lw = 0.4)  
  boxplot (aggr.region$Coral_cover ~ aggr.region$Year,  
           ylim = c(0, 100),           
           axes = FALSE,
           cex = 0.5,                 
           lwd = 1,
           pars = list(boxwex = 0.7, staplewex = 0.3, outwex = 0.3),
           las = 2,                  
           cex.axis=0.5,
           xlab = "Year",
           ylab = "") 
  axis(2,lwd = 0.5)
  box(lwd = 0.5)
  x <-c(1:45)
  YearsR = as.character(c(1970:2014))
  axis(1,  at = x, labels = YearsR, lwd = 0.5, cex.axis=0.5, las = 2)

  dev.off()

## null device 
##           1

Polynomial adjustment for all regions

library(ggplot2)

theme_set (theme_bw() + theme(panel.border = element_blank(), 
                              panel.grid.major = element_blank(),
                              panel.grid.minor = element_blank(), 
                              axis.line = element_line(colour = "black")))
x = aggr.region$Year
y = aggr.region$Coral_cover
qplot(x,y, geom='smooth', span = 0.3)+ xlim(1970, 2014)+xlab("Year")+
   ylab("Coral cover")+ylim(0,80)+ geom_point()

## Warning: Removed 1 rows containing non-finite values (stat_smooth).

## Warning: Removed 1 rows containing missing values (geom_point).

DHW and coral cover change

# Load packages 
library(tidyverse)
library(dplyr)
library(data.table)
library(RColorBrewer)

FIGURE 2A. The sequence of yearly maximum DHW experienced by the best-represented coral reefs sites from 1982 to 2014

Select and transform DHW data

# Select DHW data file: maxDHW.csv
  dhw <- read.csv("Data/maxDHW.csv", header=TRUE)

# Exclude years 1981, and 2017-2019
  dhw <- dhw[!dhw$Year == 1981, ] 
  dhw <- filter(dhw, Year < 2015)

# Create groups for plot
  ts.Cano_island <- subset(dhw, dhw$Site == "Cano_island")     
  ts.Uva_island <- subset(dhw, dhw$Site == "Uva_Island")          
  ts.Gorgona_island <- subset(dhw, dhw$Site == "Gorgona_Island")          
  ts.Darwin_island <- subset(dhw, dhw$Site == "Darwin_North_Anchorage")          
  ts.Bahia_Banderas <- subset(dhw, dhw$Site == "Bahia_Banderas")

Plot max DHW in each location

# Year, maxDHW
  par(oma = c(1,1,1,1)) 
  par(mar = c(4,4,1,1), cex.axis = 1,  las = 1)  
  par(lwd = 1)   
  
  plot(ts.Cano_island$Year,ts.Cano_island$maxDHW,   # Create a white false line with real data
       xlab="Years", ylab="Degree heating weeks (°C-weeks)",  cex=1.5, cex.lab=0.6, col="white",axes = FALSE,
       ylim = c(0, 28), lwd = 1.5) 
  par(new=T)
  rect(1982, 0, 1983.5, 28, col="#fee8c8", border ="#fee8c8")           # Strong ENSO years
  par(new=T)
  rect(1997, 0, 1998.5, 28, col="#fee8c8", border ="#fee8c8")
  par(new=T)
  plot(ts.Cano_island$Year, ts.Cano_island$maxDHW,type="l",              # Thermally stable
  xlab="", ylab="", lwd = 1.5, col="#8c2d04",axes = FALSE,
       ylim = c(0, 28))
  par(new=T)
  plot(ts.Uva_island$Year, ts.Uva_island$maxDHW, type="l",                # Thermally stable
       xlab="", ylab="", lwd = 1.5, col="#feb24c",axes = FALSE,
       ylim = c(0, 28))
  par(new=T)
  plot(ts.Gorgona_island$Year,ts.Gorgona_island$maxDHW, type="l",         # Tropical upwelling
       xlab="", ylab="",lwd = 1.5, col="#fb6a4a", axes = FALSE,
       ylim = c(0, 28))
  par(new=T)
  plot(ts.Darwin_island$Year,ts.Darwin_island$maxDHW, type="l",           # Galapagos
       xlab="", ylab="", lwd = 1.5, col="#9e9ac8",axes = FALSE,
       ylim = c(0, 28))
  par(new=T)
  plot(ts.Bahia_Banderas$Year, ts.Bahia_Banderas$maxDHW, type="l",        # Seasonal
       xlab="", ylab="",lwd = 1.5,col="#a6d96a", axes = FALSE,
       ylim = c(0, 28))
  y <- c(0,4,8,12,16,20,24,28)                                            # Axis
  axis(2,lwd = 0.5, at = y,cex.axis=0.6)
  box(lwd = 0.5)
  x <-c(1982:2014)
  YearsR = as.character(c(1982:2014))
  axis(1,  at = x, labels = YearsR, lwd = 0.4, cex.axis=0.5, las = 2) 
  abline(h=4, col="#525252", lty=2)
  abline(h=8, col="#525252", lty=2)

# Legend
  sites <- c("Thermally stable - Caño Island","Thermally stable - Uva Island",
             "Tropical upwelling - Gorgona Island","Equatorial Upwelling - Darwin Island", 
             "Seasonal - Bahía Banderas")
  par(xpd=TRUE)     
  legend(x = 2000, y = 20, 
         xpd = NA,  # or TRUE 
         inset=c(-0.5,0),
         cex = 0.4,
         legend = sites,         # box categories 
         col = c("#8c2d04","#feb24c","#fb6a4a","#9e9ac8","#a6d96a"),
         lty= 1,
         lwd = 1.5,
         x.intersp = 0.5,
         y.intersp = 2, 
         bty="n")

  dev.off() # uncomment to reset parameters for plot

## null device 
##           1

Relationship of coral cover rate of change and thermal stress

The DHW index measures the accumulated thermal stress above a bleaching threshold for 12 consecutive weeks.
At each site, the maximum monthly mean sea surface # temperature (MMM) is calculated from a long-term satellite climatology dataset.
Any temperature surpassing the MMM by 1ºC or more is added for 12 consecutive weeks,
because temperatures ~1ºC above the MMM are known to cause coral bleaching (Liu et al., 2006).
For example, if the water temperature during four consecutive weeks exceeds the MMM in 1.3°C, # 1.1°C, 1.2°C, and 1.4°C, then,
the site accumulates 5°C-weeks when these anomalies are summed.

We tested the hypothesis that the strength of the Log_annual_rate_of_change (response variable) is a function of the maxDHW experienced. We had measured the rate of change from four Thermal regimes, and fitted Thermal regimes as a random intercept, and estimate a common slope (change in the rate of change) for maxDHW across all Thermal regimes by fitting it as a fixed effect.

## Load data
  roc.cover <- read.csv("Data/RoC_DHW.csv", header=TRUE)  ### Select file:  RoC_DWH.csv

# Load packages
  library(lme4)
  library(lmerTest)
  library(MuMIn)

Model 6: coral cover change and DHW

model6 <- lmerTest::lmer(Log_annual_rate_of_change ~ maxDHW + (1|Thermal_regime), data=roc.cover)
step(model6)

## Backward reduced random-effect table:
## 
##                      Eliminated npar  logLik    AIC    LRT Df Pr(>Chisq)
## <none>                             4 -17.738 43.476                     
## (1 | Thermal_regime)          0    3 -27.772 61.545 20.069  1  7.469e-06
##                         
## <none>                  
## (1 | Thermal_regime) ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##        Eliminated  Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## maxDHW          0 0.65054 0.65054     1 128.59   9.873 0.002082 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## Log_annual_rate_of_change ~ maxDHW + (1 | Thermal_regime)

summary(model6)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Log_annual_rate_of_change ~ maxDHW + (1 | Thermal_regime)
##    Data: roc.cover
## 
## REML criterion at convergence: 35.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.7962 -0.1472  0.1167  0.4152  1.8621 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  Thermal_regime (Intercept) 0.02501  0.1581  
##  Residual                   0.06589  0.2567  
## Number of obs: 133, groups:  Thermal_regime, 4
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)   
## (Intercept)  -0.085763   0.084256   3.155381  -1.018  0.38032   
## maxDHW       -0.016316   0.005193 128.594887  -3.142  0.00208 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##        (Intr)
## maxDHW -0.165

r.squaredGLMM(model6)   #  returns the marginal and the conditional R2

## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help
## page.

##             R2m       R2c
## [1,] 0.05249512 0.3131505

# marginal R2  describes the proportion of variance explained by the fixed factor(s) alone:
# conditional R2 describes the proportion of variance explained by both the fixed and random factors
# Nakagawa, S. & Schielzeth, H. (2013).

FIGURE 2B: Relationship of coral cover rate of change and thermal stress

# Load packages
  library(tidyverse)
  library(dplyr)
  library(data.table)
  library(RColorBrewer)
 
# Create groups for plot
  Seasonal          = subset(roc.cover, roc.cover$Thermal_regime == "Seasonal")          
  Thermally_Stable  = subset(roc.cover, roc.cover$Thermal_regime == "Thermally_Stable")  
  Tropical_Upwelling = subset(roc.cover, roc.cover$Thermal_regime == "Upwelling") 
  Galapagos         = subset(roc.cover, roc.cover$Thermal_regime == "Galapagos_Islands") 

# Plot
  par(oma = c(1,1,1,1)) 
  par(mar = c(4,4,1,1), cex.axis = 1,  las = 1) 
  par(lwd = 1)
  plot (Seasonal$maxDHW, Seasonal$Log_annual_rate_of_change,
        pch = 21, lwd = 0.5, bg = "#a6d96a93", col="black", axes = FALSE,
        cex=1.5,cex.lab=0.6,
        xlim = c(0,30), ylim = c(-1.75, 1),
        xlab = "Degree heating weeks (°C-weeks)",
        ylab = "Annual rate of change in coral cover")
  par(new=T); 
  plot (Thermally_Stable$maxDHW, Thermally_Stable$Log_annual_rate_of_change,
        xlab="", ylab="", pch = 21, lwd = 0.5, col="black", bg = "#feb24c93", axes = FALSE,
        cex=1.5,cex.lab=0.6,
        xlim = c(0,30), ylim = c(-1.75, 1))
  par(new=T); 
  plot (Tropical_Upwelling$maxDHW, Tropical_Upwelling$Log_annual_rate_of_change,
        xlab="", ylab="", pch = 21, lwd = 0.5, bg = "#fb6a4a93",col="black", axes = FALSE,
        cex=1.5,cex.lab=0.6,
        xlim = c(0,30), ylim = c(-1.75, 1))
  par(new=T); 
  plot (Galapagos$maxDHW, Galapagos$Log_annual_rate_of_change,
        xlab="", ylab="", pch = 21, lwd = 0.5, bg = "#9e9ac893",col="black", axes = FALSE,
        cex=1.5,cex.lab=0.6,
        xlim = c(0,30), ylim = c(-1.75, 1))
  box(lwd = 0.5)
  x <-c(0,5,10,15,20,25,30)
  axis(1, lwd = 0.5, at = x, cex.axis = 0.5)
  y <- c(-1.5,-1.0,-0.5,0,0.5,1)
  axis(2,lwd = 0.5, at = y, cex.axis = 0.5)
  abline(v=4, col="#525252", lty=2); abline(v=8, col="#525252", lty=2); abline(h=0, col="#525252", lty=2); 
  # Axis text I, II, III, IV
  text(3.5, 0.3, "II", cex = 0.6); text(5, 0.3, "I", cex = 0.6)
  text(3.5, - 0.3, "III", cex = 0.6); text(5, -0.3, "IV", cex = 0.6)
  
  par(xpd=TRUE) 
  legend(x = 10, y = -0.5, 
         xpd = NA,   
         inset=c(-0.5,0),
         legend = c("Equatorial Upwelling", "Thermally stable", "Tropical upwelling","Seasonal"),
         pch=21,
         pt.bg = c("#9e9ac8","#fb6a4a","#fb6a4a","#a6d96a"),
         x.intersp = 0.5,
         y.intersp = 2, 
         bty="n")

  dev.off() # uncomment to reset parameters for plot

## null device 
##           1

Descriptives statistics heat stress

No stress

no.stress <- subset(roc.cover, maxDHW < 4)                      #subset of DHW < 4
length(no.stress$maxDHW)/length(roc.cover$maxDHW)               # ratio no stress: total

## [1] 0.8496241

1- length(no.stress$maxDHW)/length(roc.cover$maxDHW)               # ratio stress: total

## [1] 0.1503759

positive.rate <- subset(no.stress, Log_annual_rate_of_change > 0)   #subset of Rate > 0
length(positive.rate$Log_annual_rate_of_change)                     # Number of observations in

## [1] 48

negative.rate <- subset(no.stress, Log_annual_rate_of_change < 0)   #subset of Rate > 0
length(negative.rate$Log_annual_rate_of_change)                     # Number of observations in

## [1] 58

No stress seasonal

no.stress.seasonal <- subset(no.stress, Thermal_regime == "Seasonal")   
positive.rate <- subset(no.stress.seasonal, Log_annual_rate_of_change > 0)  
negative.rate <- subset(no.stress.seasonal, Log_annual_rate_of_change < 0)  
length(positive.rate$Log_annual_rate_of_change)/length(negative.rate$Log_annual_rate_of_change)

## [1] 0.07142857

No stress upwelling

no.stress.upwelling <- subset(no.stress, Thermal_regime == "Upwelling")
positive.rate <- subset(no.stress.upwelling, Log_annual_rate_of_change > 0)  
negative.rate <- subset(no.stress.upwelling, Log_annual_rate_of_change < 0)  
length(positive.rate$Log_annual_rate_of_change)/length(negative.rate$Log_annual_rate_of_change)

## [1] 1

No stress thermally stable

no.stress.Thermally_Stable <- subset(no.stress, Thermal_regime == "Thermally_Stable")   
positive.rate <- subset(no.stress.Thermally_Stable, Log_annual_rate_of_change > 0)  
negative.rate <- subset(no.stress.Thermally_Stable, Log_annual_rate_of_change < 0)  
length(positive.rate$Log_annual_rate_of_change)/length(negative.rate$Log_annual_rate_of_change)

## [1] 1.2

No stress Galapagos

no.stress.Galapagos_Islands <- subset(no.stress, Thermal_regime == "Galapagos_Islands") 
positive.rate <- subset(no.stress.Galapagos_Islands, Log_annual_rate_of_change > 0)  
negative.rate <- subset(no.stress.Galapagos_Islands, Log_annual_rate_of_change < 0)  
length(positive.rate$Log_annual_rate_of_change)/length(negative.rate$Log_annual_rate_of_change)

## [1] 1

Supplementary figures

library(pals)
library(dplyr)
library(reshape2)
library(RColorBrewer)

FIGURE S2: a) Number of regions surveyed per year, b): Number of countries surveyed per year

# dev.off() # uncomment to reset parameters for plot

# Barplot number of surveys per year
  obs_num = as.data.frame(table(cover$Year,cover$Region))    
  colnames(obs_num) = c("Year", "Region", "Freq")
  wide_obs_num = dcast(obs_num, Region~Year, value = "Freq")

## Using Freq as value column: use value.var to override.

  df = data.matrix(wide_obs_num[,2:45])   
  as.table(df)

##   1970 1971 1972 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984
## A    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## B    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## C    0    0    1    0    0    0    0    0    1    1    0    0    1    0
## D    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## E    0    0    0    0    0    0    0    0    0    1    0    1    0    4
## F    0    0    0    2    3    2    0    0    0    0    0    1    0    1
## G    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## H    3   10    2   13    4    4    3    5    2    5    2    3    5    3
## I    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## J    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## K    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## L    0    0    0    0    0    0    0    2    0    0    0    0    0    0
## M    0    0    0    0    0    0    0    0    0    0    0    1    0    0
## N    2    6    0    2    1    0    0    0    0    0    0    1    4    2
## O    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## P    0    0    0    0    0    0    0    0    0    0    0    0    0    0
##   1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998
## A    0    0    0    0    0    0    0    0    0    1    0    0    0    0
## B    0    0    5    0    0    0    0    1    0    0    0    0    0    0
## C    0    0    1    2    4    0    0    1    0    0    2    2    0    4
## D    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## E   12    1    2    0    1    0    0   11    0    2    0    6    0    1
## F    0    1    1    0    1    0    0    0    2    0    0    0    0    0
## G    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## H    4    5    2    3    1    1    0    4    1    4    4    0    1    1
## I    0    0    0    0    0    0    2    0    0    0    0    0   16   10
## J    0    0    1    1    1    0    3    1    1    0    0    0    2    2
## K    0    0    0    0    0    0    0    0    0    0    0    0    2    1
## L    0    0    0    0    0    0    0    0    0    1    1    4    4    2
## M    0    0    0    0    0    0    0    1    0    0    0    0    0    0
## N    0    1    5    1    2    0    1    1    1    1    1    1    1    1
## O    0    0    0    0    0    0    0    0    1    1    0    0    0    0
## P    0    0    0    0    0    0    0    0    0    0    0    0    0    0
##   1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012
## A    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## B    0    0    0    1    0    0    0    1    0    0    4    0    0    0
## C    7    4    3    4    6    4    3    6    3    2    3    2    2    3
## D    0    0    0    0    0    0    0    0    0    0    0    4    0    1
## E    1    0    2    0    4    2    2    0    1    0    2    0    0    0
## F    0    0    1    1    0    1    3    6    2    2    0    0    1    0
## G    0    0    0    0    0    0    0    0    0    0    1    0    0    0
## H    0    1    4   14   12    1    1    1    1    0    1    2    0    0
## I    0    0    0    6    3    0    1    1    4    0   17    6    2    0
## J    0    0    0    1    2    0    0    0    0    5    7    0    0    1
## K    1    5    5    3    2    0    0    0    0    0    2    0    0    0
## L    0    0    2    0    1    0    1    3    2    0   17    3    1    1
## M    0    0    0    0    0    1    2   14    8    0    0    0    0    1
## N    1    1    1    1    1    1    0    0    0    0    1    0    0    0
## O    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## P    0    0    0    0    0    2    0    0    0    0    0    0    0    0
##   2013 2014
## A    0    0
## B    1    1
## C    2    1
## D    4    2
## E    0    0
## F    0    0
## G    0    0
## H    0    0
## I    2    2
## J    0    0
## K    0    0
## L    1    3
## M    0    0
## N    0    0
## O    0    0
## P    0    0

  df[ df>1 ] <- 1 
  levels(cover$Region) # How many regions? = how many colours

##  [1] "Clipperton"                 "Cocos_Islands"             
##  [3] "Colombia"                   "Costa_Rica_Central"        
##  [5] "Costa_Rica_Southern"        "Eastern_Galapagos_Islands" 
##  [7] "El_Salvador"                "Gulf_of_Chiriqui"          
##  [9] "Mexican_Central"            "Mexican_Northern"          
## [11] "Mexican_Southern"           "Nicaragua_Papagayo_Zone"   
## [13] "Northern_Galapagos_Islands" "Panama_Bight"              
## [15] "Revillagigedo_Islands"      "Western_Galapagos_Islands"

  Regions <- c(levels(obs_num$Region))
  color.regions = warmcool(17)
  
  # Figure 2a
  par(mfcol=c(2,1), cex.main = 1.5, mar = c(4,4,1,1), cex.lab = 1.5, cex.axis = 1,  las = 1)
  barplot(df,
          ylim=c(0,15),
          col = color.regions,   
          ylab="Number of regions surveyed per year",  
          xlab="Year",
          #lwd = 1,          
          las = 2,           
          axis.lty = 1, 
          cex.names = 0.4,   
          cex.axis = 0.8)
# Legend as a box
  legend( #"bottomright",    
     x = 2, y = 16, 
     xpd = NA,  # or TRUE 
     legend = Regions,          
     fill = color.regions[1:17],   
     inset=c(-0.5,0),
     cex = 0.35,
     text.font=0.5,            
     x.intersp = -0.5,
     y.intersp = 1,
     lwd = 0.1,
     bty="n")

# Figure 2b
obs_num = as.data.frame(table(cover$Year,cover$Country))    # Table per country
colnames(obs_num) = c("Year", "Country", "Freq")
wide_obs_num = dcast(obs_num, Country~Year, value = "Freq")

## Using Freq as value column: use value.var to override.

dfc = data.matrix(wide_obs_num[,2:45])   
as.table(dfc)

##   1970 1971 1972 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984
## A    0    0    1    0    0    0    0    0    1    1    0    0    1    0
## B    0    0    2    9    0    0    0    2    0    2    0    1    0    4
## C    0    0    0    2    3    2    0    0    0    0    0    2    0    1
## D    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## E    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## F    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## G    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## H    5   16    0    6    5    4    3    5    2    4    2    4    9    5
##   1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998
## A    0    0    1    2    4    0    0    1    0    0    2    2    0    4
## B   12    1    7    0    1    0    0   15    0    6    4   10    4    3
## C    0    1    1    0    1    0    0    1    2    0    0    0    0    0
## D    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## E    0    0    0    0    0    0    0    0    0    1    0    0    0    0
## F    0    0    1    1    1    0    5    1    2    1    0    0   20   13
## G    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## H    4    6    7    4    3    1    1    2    2    2    2    1    2    2
##   1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012
## A    7    4    3    4    6    4    3    6    3    2    3    2    2    3
## B    1    0    7    1    5    2    3    4    3    0   12    7    1    2
## C    0    0    1    1    0    4    5   20   10    2    0    0    1    1
## D    0    0    0    0    0    0    0    0    0    0    1    0    0    0
## E    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## F    1    5    5   10    7    0    1    1    4    5   26    6    2    1
## G    0    0    0    0    0    0    0    0    0    0   11    0    0    0
## H    1    2    2   15   13    2    1    1    1    0    2    2    0    0
##   2013 2014
## A    2    1
## B    6    6
## C    0    0
## D    0    0
## E    0    0
## F    2    2
## G    0    0
## H    0    0

dfc[ dfc>1 ] <- 1 # replace all values greater than one with one
color.countries <- brewer.pal(9, "Oranges")
countries = c(levels(obs_num$Country))
barplot(dfc,
        col = color.countries,   
        ylim=c(0,10),
        ylab="Number of countries surveyed per year",  
        xlab="Year",
        las = 2,           
        axis.lty = 1,      
        cex.names = 0.4,   
        cex.axis = 0.8)
# Legend as a box
legend( #"bottomright",    
   x = 2, y = 10, 
   xpd = NA,  # or TRUE 
   legend = countries,         
   fill = color.countries[1:9],   
   inset=c(-0.5,0),
   cex = 0.35,
   text.font=1,            
   x.intersp = -0.5,
   y.intersp = 1,
   lwd = 0.1,
   bty="n")

FIGURE S3: Variation in the live coral cover for 1970–2014 per thermal regime

Temporal variation in the live coral cover for 1970–2014. * (a) Tropical upwelling regime. * (b) Seasonal regime. * (c) Thermally stable regime. * (d) Equatorial upwelling regime. The smooth line represents the best fit to the time series with a 95% confidence level interval and a span = 0.3.

# dev.off() # uncomment to reset parameters for plot

library(ggplot2)
library(lemon)

aggr.region$Thermal_regime<-factor(aggr.region$Thermal_regime, 
                                      levels=c("Upwelling", "Seasonal", "Thermally_Stable", "Galapagos_Islands"), 
                                      labels=c("(a)", "(b)", "(c)", "(d)"))
summary(aggr.region$Thermal_regime)

## (a) (b) (c) (d) 
##  75  28  74  25

ggplot(aggr.region, aes(x=Year, y=Coral_cover)) + 
  theme_bw() + theme(panel.border = element_blank(),   # set ggplot to white background
                              panel.grid.major = element_blank(),
                              panel.grid.minor = element_blank(), 
                              axis.line = element_line(),
                              strip.background = element_blank(),
                              strip.text = element_text(face="bold",hjust = 0, vjust = -1))+
  geom_point()+ geom_smooth (span = 0.3) +
  scale_x_continuous(limits=c(1970, 2014), expand = c(0.01, 0.01), "Year") +
  scale_y_continuous(limits=c(-19, 91), expand = c(0, 0), "Coral cover") +
  facet_wrap(Thermal_regime~., ncol=1, scales = "free_x", strip.position ="top")

Coral cover analysis

Mauricio Romero-Torres

September 8, 2019